import numpy as np
import json
import math
import gurobipy as gp
from gurobipy import GRB
import pandas as pd
from datetime import datetime, timedelta


def CreateEvConstraint(SAModel, periods, n_EV, Ev_n, EVPmax,EVsChargingStartTime, EVsChargingEndTime, EnergyIni,EVsIniChargePower,Mode):
    EVsOptChargePower = SAModel.addVars(n_EV, periods, ub=EVPmax, vtype=GRB.CONTINUOUS)  # 各个电动汽车各时段充电功率变量
    '''未接网时段功率为0，接网时段充电量在一定限度内'''
    for i in range(n_EV):
        if EVsChargingStartTime[i] < EVsChargingEndTime[i]:
            for T in range(EVsChargingStartTime[i]):
                SAModel.addConstr(EVsOptChargePower[i, T] == 0)
            for T in range(EVsChargingEndTime[i], periods):
                SAModel.addConstr(EVsOptChargePower[i, T] == 0)
            # 计算这段时间的充电电量
            EnergyOpt = gp.quicksum(
                EVsOptChargePower[i, T] * Ev_n for T in range(EVsChargingStartTime[i], EVsChargingEndTime[i])) / (
                                   periods / 24)
            SAModel.addConstr(EnergyOpt == EnergyIni[i])
        else:
            for T in range(EVsChargingEndTime[i], EVsChargingStartTime[i]):
                SAModel.addConstr(EVsOptChargePower[i, T] == 0)
            # 计算这段时间的充电电量
            EnergyOpt = gp.quicksum(EVsOptChargePower[i, T] * Ev_n for T in range(EVsChargingStartTime[i], periods)) / (
                        periods / 24)  \
                       + gp.quicksum(EVsOptChargePower[i, T] * Ev_n for T in range(EVsChargingEndTime[i])) / (
                                   periods / 24)
            SAModel.addConstr(EnergyOpt == EnergyIni[i])
        if Mode == 2:
            '''只调AC，EV不动，保持原功率'''
            for T in range(periods):
                SAModel.addConstr(EVsOptChargePower[i, T] == EVsIniChargePower[i, T])


    return EVsOptChargePower

def CreateAcConstraint(SAModel, periods, n_AC, EnergyIni,Ca, Ra, Cm, Rm, To, Ta0, Tm0, Tmax, Tmin, Cope,ACsIniChargePower,Mode):
    Pac = SAModel.addVars(n_AC, periods, lb=-GRB.INFINITY, ub=0, vtype=GRB.CONTINUOUS)
    for i in range(n_AC):
        SAModel.addConstr(-gp.quicksum(Pac[i,k] for k in range(periods)) == EnergyIni[i])

    Ta = SAModel.addVars(n_AC, periods, lb=-GRB.INFINITY, vtype=GRB.CONTINUOUS)
    Tm = SAModel.addVars(n_AC, periods, lb=-GRB.INFINITY, vtype=GRB.CONTINUOUS)
    step = 24 / periods
    for i in range(n_AC):
        # 计算解析递推形式的参数
        b = 1 / (Ca[i] * Ra[i]) + 1 / (Ca[i] * Rm[i]) + 1 / (Cm[i] * Rm[i])
        c = 1 / (Ca[i] * Cm[i] * Ra[i] * Rm[i])
        delta = b ** 2 - 4 * c
        if delta > 0:
            A = np.zeros((periods, 2, 2))
            B = np.zeros((periods, 2))
            r1 = (-b + math.sqrt(b ** 2 - 4 * c)) / 2
            r2 = (-b - math.sqrt(b ** 2 - 4 * c)) / 2

            e1 = math.exp(r1 * step)
            e2 = math.exp(r2 * step)
            for k, t in enumerate(np.arange(0, 24, step)):
                A[k, 0, 0] = (r1 * e1 - r2 * e2) / math.sqrt(delta) + 1 / (Cm[i] * Rm[i]) * (e1 - e2) / math.sqrt(delta)
                A[k, 0, 1] = 1 / (Ca[i] * Rm[i]) * (e1 - e2) / math.sqrt(delta)
                A[k, 1, 0] = 1 / (Cm[i] * Rm[i]) * (e1 - e2) / math.sqrt(delta)
                A[k, 1, 1] = (r1 * e2 - r2 * e1) / math.sqrt(delta) - 1 / (Cm[i] * Rm[i]) * (e1 - e2) / math.sqrt(delta)
                B[k, 0] = 1 + ((r2 * e1 - r1 * e2) + 1 / (Ca[i] * Ra[i]) * (e1 - e2)) / math.sqrt(delta)
                B[k, 1] = 1 + (r2 * e1 - r1 * e2) / math.sqrt(delta)

        else:
            A = np.zeros((periods, 2, 2))
            B = np.zeros(periods, 2)
            r = -b / 2
            w = math.sqrt(-delta) / 2
            for k, t in enumerate(np.arange(0, 24, step)):
                A[k, 0, 0] = math.exp(r * t) * math.cos(w * t) + math.exp(r * t) * math.sin(w * t) / w * (
                        1 / (Cm[i] * Rm[i]) + r)
                A[k, 0, 1] = math.exp(r * t) * math.sin(w * t) / w / (Ca[i] * Rm[i])
                A[k, 1, 0] = math.exp(r * t) * math.sin(w * t) / w / (Cm[i] * Rm[i])
                A[k, 1, 1] = math.exp(r * t) * math.cos(w * t) - math.exp(r * t) * math.sin(w * t) / w * (
                        1 / (Cm[i] * Rm[i]) + r)

                B[k, 0] = 1 - math.exp(r * t) * math.cos(w * t) + math.exp(r * t) * math.sin(w * t) / w * (
                        1 / (Ca[i] * Ra[i]) + r)
                B[k, 1] = 1 - math.exp(r * t) * math.cos(w * t) + math.exp(r * t) * math.sin(w * t) / w * r

        # 建立温度之间的关系
        for k in range(periods):
            if k == 0:
                SAModel.addConstr(Ta[i, k] == A[k, 0, 0] * Ta0[i] + A[k, 0, 1] * Tm0[i] + B[k, 0] * (
                            To[k] + Pac[i, k] * Cope[i] * Ra[i]))
                SAModel.addConstr(Tm[i, k] == A[k, 1, 0] * Ta0[i] + A[k, 1, 1] * Tm0[i] + B[k, 1] * (
                            To[k] + Pac[i, k] * Cope[i] * Ra[i]))
            else:
                SAModel.addConstr(Ta[i, k] == A[k, 0, 0] * Ta[i, k - 1] + A[k, 0, 1] * Tm[i, k - 1] + B[k, 0] * (
                            To[k] + Pac[i, k] * Cope[i] * Ra[i]))
                SAModel.addConstr(Tm[i, k] == A[k, 1, 0] * Ta[i, k - 1] + A[k, 1, 1] * Tm[i, k - 1] + B[k, 1] * (
                            To[k] + Pac[i, k] * Cope[i] * Ra[i]))

            SAModel.addConstr(Ta[i, k] <= Tmax)
            SAModel.addConstr(Ta[i, k] >= Tmin)
            # SAModel.addConstr(Pac[i,k] <= Pmax)
        if Mode == 3:
            '''只调EV，AC不动，保持原功率'''
            for T in range(periods):
                SAModel.addConstr(-Pac[i, T] == ACsIniChargePower[i, T])
    return Pac


'''
breif:
    优化函数
Params:
    DERData:从json文件中读取出的资源参数
    IniData:从数据库中读取得到的Dataframe数据
Return:
    各资源优化后的功率序列数据（Dataframe类型）
'''
def VceOpt(DERData,IniData):
    Mode = DERData['Mode']
    Tmax = DERData['AC_Tmax']
    Tmin = DERData['AC_Tmin']
    AreasId = list(DERData.keys())[3:]
    AreasOptPower = {}  # 存放各个台区的优化变量
    MultiAreaModel = gp.Model(name='Multi')
    periods = 24
    for n, area in enumerate(AreasId):
        AreasOptPower[area] = {}
        Data = DERData[area] # 读取当前台区的数据
        '''整合空调参数'''
        n_AC = len(Data['ACsParam'])
        To = np.array(IniData[f'{area}_Temperature'])  # 温度属于时间序列数据，从dataframe中读取
        Ca = np.array([Data['ACsParam'][f'AC_{i}']['Ca'] for i in range(n_AC)])
        Ra = np.array([Data['ACsParam'][f'AC_{i}']['Ra'] for i in range(n_AC)])
        Cm = np.array([Data['ACsParam'][f'AC_{i}']['Cm'] for i in range(n_AC)])
        Rm = np.array([Data['ACsParam'][f'AC_{i}']['Rm'] for i in range(n_AC)])
        Cope = np.array([Data['ACsParam'][f'AC_{i}']['Cope'] for i in range(n_AC)])
        Ta0 = np.array([Data['ACsParam'][f'AC_{i}']['Ta0'] for i in range(n_AC)])
        Tm0 = np.array([Data['ACsParam'][f'AC_{i}']['Tm0'] for i in range(n_AC)])
        ACsEnergyIni = np.array([Data['ACsParam'][f'AC_{i}']['EnergyIni'] for i in range(n_AC)])
        '''整合EV参数'''
        n_EV = len(Data['EVParam'])
        Ev_n = Data['EVParam']['EV_0']['Ev_n']
        EVPmax = Data['EVParam']['EV_0']['ChargePowerMax']
        EVsChargingStartTime = np.array([Data['EVParam'][f'EV_{i}']['ChargingStartTime'] for i in range(n_EV)])
        EVsChargingEndTime = np.array([Data['EVParam'][f'EV_{i}']['ChargingEndTime'] for i in range(n_EV)])
        EVsSocIni = np.array([Data['EVParam'][f'EV_{i}']['EVsSocIni'] for i in range(n_EV)])
        EVsEnergyIni = np.array([Data['EVParam'][f'EV_{i}']['EnergyIni'] for i in range(n_EV)])

        AreasOptPower[area]['n_EV'] = n_EV
        AreasOptPower[area]['n_AC'] = n_AC
        AreasOptPower[area]['To'] = To
        '''从Inidata里读取EV和AC的原始功率曲线'''
        EVsIniChargePower = np.zeros((n_EV,periods))
        ACsIniChargePower = np.zeros((n_AC, periods))
        for i in range(n_EV):
            EVsIniChargePower[i] = np.array(IniData[f'{area}_EV_{i}_Ini'])
        for i in range(n_AC):
            ACsIniChargePower[i] = np.array(IniData[f'{area}_AC_{i}_Ini'])

        AreasOptPower[area]['EVsOptChargePower'] \
            = CreateEvConstraint(MultiAreaModel, periods, n_EV, Ev_n, EVPmax,
                                     EVsChargingStartTime, EVsChargingEndTime, EVsEnergyIni, EVsIniChargePower, Mode)

        AreasOptPower[area]['ACsOptChargePower'] \
            = CreateAcConstraint(MultiAreaModel, periods, n_AC, ACsEnergyIni, Ca, Ra, Cm, Rm, To, Ta0, Tm0, Tmax, Tmin,
                                     Cope, ACsIniChargePower, Mode)

        AreasOptPower[area]['PVPower'] = np.array(IniData[f'{area}_PV'])  # 光伏属于时间序列数据，从dataframe中读取
        AreasOptPower[area]['LDPower'] = np.array(IniData[f'{area}_LD'])  # 负荷属于时间序列数据，从dataframe中读取

        AreasOptPower[area]['NetPower'] = MultiAreaModel.addVars(periods, lb=-GRB.INFINITY, vtype=GRB.CONTINUOUS)
        AreasOptPower[area]['Pmax'] = MultiAreaModel.addVar(vtype=GRB.CONTINUOUS)
        for T in range(periods):
            MultiAreaModel.addConstr(AreasOptPower[area]['NetPower'][T] ==
                                     gp.quicksum(
                                         AreasOptPower[area]['EVsOptChargePower'][i, T] for i in range(n_EV)) + gp.quicksum(
                -AreasOptPower[area]['ACsOptChargePower'][i, T] for i in range(n_AC)) + AreasOptPower[area]['LDPower'][T] -
                                     AreasOptPower[area]['PVPower'][T])
            MultiAreaModel.addConstr(AreasOptPower[area]['Pmax'] >= AreasOptPower[area]['NetPower'][T])
            MultiAreaModel.addConstr(AreasOptPower[area]['Pmax'] >= -AreasOptPower[area]['NetPower'][T])

    Pmain = MultiAreaModel.addVars(periods, lb=-GRB.INFINITY, vtype=GRB.CONTINUOUS)
    Pmain_max = MultiAreaModel.addVar(lb=0, vtype=GRB.CONTINUOUS)
    for T in range(periods):
        MultiAreaModel.addConstr(Pmain[T] == sum(AreasOptPower[area]['NetPower'][T] for n, area in enumerate(AreasId)))
        MultiAreaModel.addConstr(Pmain_max >= Pmain[T])
        MultiAreaModel.addConstr(Pmain_max >= -Pmain[T])

    obj = Pmain_max + sum(AreasOptPower[area]['Pmax'] for n, area in enumerate(AreasId))
    MultiAreaModel.setObjective(obj, GRB.MINIMIZE)
    MultiAreaModel.optimize()
    '''优化结束取出各个资源的功率优化结果'''
    DERsPowerList = []  # 存放各个台区的资源功率序列，并最终将其合并
    PowerRowIndexList = []
    for m, area in enumerate(AreasId):
        ACPower = np.array([[-AreasOptPower[area]['ACsOptChargePower'][i, T].x for T in range(periods)] for i in
                            range(AreasOptPower[area]['n_AC'])])
        EVPower = np.array([[AreasOptPower[area]['EVsOptChargePower'][i, T].x for T in range(periods)] for i in
                            range(AreasOptPower[area]['n_EV'])])
        PVPower = AreasOptPower[area]['PVPower']
        LDPower = AreasOptPower[area]['LDPower']
        To = AreasOptPower[area]['To']
        Power = np.vstack((ACPower, EVPower))
        Power = np.vstack((Power, PVPower))
        Power = np.vstack((Power, LDPower))
        Power = np.vstack((Power, To))
        DERsPowerList.append(Power)
        PowerRowIndexList.append([area + '_' + f'AC_{i}' + '_Opt' for i in range(AreasOptPower[area]['n_AC'])]
                                 + [area + '_' + f'EV_{i}' + '_Opt' for i in range(AreasOptPower[area]['n_EV'])]
                                 + [area + '_' + 'PV'] + [area + '_' + 'LD'] + [area + '_' + 'Temperature'])
    DERsPower = np.vstack(DERsPowerList)  # 按顺序存放每一个资源的初始功率曲线
    PowerRowIndex = sum(PowerRowIndexList, [])  # 生成行名称，即每一个资源的名称
    # 生成时间列名称，从2024-10-18 00:00开始，间隔1小时
    start_time = datetime(2024, 10, 18, 0, 0)
    time_intervals = [start_time + timedelta(hours=i) for i in range(24)]
    TimeColumnIndex = [t.strftime('%Y-%m-%d %H:%M') for t in time_intervals]  # 列名称为时间格式

    # 将优化结果变为一个DataFrame
    df = pd.DataFrame(DERsPower, index=PowerRowIndex, columns=TimeColumnIndex)
    df = df.T # 将其转置

    return  df

if __name__ == '__main__':
    '''---------------------------读取DataFrame类型的数据----------------------'''
    IniData = pd.read_csv('IniPower.csv')
    names = IniData.columns
    names_list = names.tolist()[1:]  # 资源名称序列
    # print(names_list)
    '''----------------------------读取柔性资源参数-----------------------------'''
    json_file_path = 'DERParams.json'
    with open(json_file_path, 'r') as json_file:
        DERData = json.load(json_file)
    '''----------------开始优化并得到优化后的各资源功率数据（Dataframe）------------'''
    OptData = VceOpt(DERData, IniData)
    
    print("功率优化结果的DataFrame创建完毕！")
    csv_file_path = 'OptPower.csv'
    OptData.to_csv(csv_file_path, index=True)  # 保存到指定路径